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We determine the form of the complex shear modulus G* in soft sphere packings near jamming. 
Viscoelastic response at finite frequency is closely tied to a packing's intrinsic relaxational modes, 
which are distinct from the vibrational modes of undamped packings. We demonstrate and explain 
the appearance of an anomalous excess of slowly relaxing modes near jamming, reflected in a di- 
verging relaxational density of states. From the density of states, we derive the dependence of G* 
on frequency and distance to the jamming transition, which is confirmed by numerics. 
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When the bubbles or droplets that make up foams 
or emulsions are packed closely together, they jam into 
a mechanically rigid state 1-5 . Near jamming they 
form amorphous packings of repulsive athermal spheres 
[2], and in recent years the linear quasistatic response 
of jammed packings has been exhaustively mapped out. 
Proximity to (un) jamming, a noncquilibrium critical 
point, organizes their elasticity, and in simulations their 
moduli scale with distance to the transition [3J. Experi- 
mental evidence for these scalings, however, remains elu- 
sive. Rheology may prove a better test bed for the jam- 
ming paradigm, as there is already a wealth of available 
data [BHS] . While there is growing numerical and theoret- 
ical evidence that jamming also organizes the rheology of 
soft sphere packings, details remain controversial [TUIU2) . 

Prior rheological studies of the jamming transition 
have focused on steady flow [HHH]- Surprisingly, de- 
spite its practical and fundamental significance, little 
is known about the oscillatory rheology of soft spheres 
near jamming. The complex shear modulus G*(uS) = 
G' (to) + iG" (ui) , composed of the storage modulus G' and 
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FIG. 1: (a) Kelvin- Voigt viscoelastic elements in series. Ele- 
ments are made of a spring with stiffness G n in parallel with a 
dashpot with coefficient r\ n \ their characterisic relaxation rate 
is s n — Gn/ri n . (b) Disk displacements of a relaxational mode 
with rate s n = 0.006 in a packing with excess coordination 
Az — 0.013; the disks are not shown, (c) Shear response for 
the same packing driven at uj = 0.006. 



loss modulus G" , describes the linear response to oscilla- 
tory driving at frequency u. It is a fundamental material 
characterization, encoding the response to shear on all 
time scales, including the quasistatic (QS) shear modu- 
lus Go in the zero frequency limit. Experiments provide 
evidence as to the form of G* : liquid and organic foams, 
emulsions, and microgel suspensions are all shear thin- 
ning, G* ~ ( tuJ ) A with A s» 0.5, over one to several 
decades in frequency pj [5] ■ 

Whereas undamped packings have vibrational normal 
modes, modes in overdamped packings are relaxational. 
In this Letter we relate relaxational modes and rheology 
to determine and explain the dependence of G* on fre- 
quency and distance to jamming. G* displays dynamical 
critical scaling; shear thinning is a critical effect and the 
shear thinning regime extends to zero frequency at the 
jamming transition. We calculate all critical exponents, 
including A = 1/2. This makes A one of the few experi- 
mentally observed jamming exponents. 

The relationship between relaxations and rheology has 
a simple schematic expression: a packing's finite fre- 
quency linear response can be mapped to a series cir- 
cuit of Kelvin- Voigt viscoelastic elements (Fig. [T^,) . Each 
element has a characteristic relaxation rate s equal to 
the eigenrate of one of the packing's relaxational modes 
(Fig. [^d). Thus the response to driving (Fig. [If;) is gov- 
erned by the distribution of relaxation rates, i.e. the den- 
sity of states D(s), which is related to, but different from, 
the vibrational density of states D(f2). Here we will de- 
termine D(s) for the first time, use scaling arguments to 
identify and explain its characteristic features, and from 
these predict the form of G*(ui). 

Model system. — Starting from numerically generated 
static packings, we study their linear response: interac- 
tions are linearized about the reference state and defor- 
mations are calculated without allowing for rearrange- 
ments of the contact network. Our jammed packings in 
d = 2 dimensions are comprised of weakly polydisperse 
disks in a L x L unit cell generated using the molecular 
dynamics protocol of Ref. [13]. Packings are character- 
ized by their mean number of contacts per particle, z. 
Static packings at the jamming transition are isostatic, 
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z = z c = 2d, hence the excess contact number Az — z—z c 
serves as a measure of distance to jamming [3l 115]. 

We impose the dynamics of the 'full bubble model' of 
Durian [TOl US] , in which foam bubbles are modeled as 
disks. Contacting disks interact via elastic and repulsive 
forces f cl = kS that are linear in the overlap S, with 
a spring constant k. Touching disks also experience a 
viscous force / Vlsc = b Av, with damping coefficient b, 
opposing their relative velocity evaluated at the contact, 
Av. The dynamics are overdamped, so that forces and 
torques balance on each bubble at every instant. We 
report stresses in units of k and times in units of b/k. 

To describe a shear deformation, which involves motion 
of the disks and distortion of the unit cell, requires 37V +1 
degrees of freedom: the 37V displacements and rotations 
of the disks, {u x ,i; u y,i> plus the magnitude 7L of the 
pure shear displacement applied to the lattice vectors of 
the unit cell. Equations of motion governing the disks' 
motions and the unit cell's deformation are most easily 
written in a Lagrangian formalism. To implement this 
we note that the elastic potential energy V is 
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where Am" and Au- 1 label normal and tangential relative 
displacements, possibly across periodic boundaries. f° 
and Ar° are the force and separation between disks in the 
reference state. Similarly, for the viscous forces described 
above, the Rayleigh dissipation function R is Q3] 
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The {pi} are disk radii. The disk rotations enter because 
relative velocities are evaluated at the contact. The forms 
of V and R will be needed below to extract the scaling 
of the relaxational density of states. 

The Lagrangian equations of motion are 



K\q(s))+sB\q(s))=a(s)L\j) 



(3) 



where the matrices are K mn = d 2 V/dq m dq n and B mn — 
d 2 R/dq m dq n , \j) is a unit vector along the strain coor- 
dinate, and a is the shear stress. We have collected all 
3N + 1 degrees of freedom in a vector \q(t)) and applied 
a Laplace transformation; transformed quantities depend 
on the independent variable s, which has units of inverse 
time and may be thought of as a relaxation rate. 

Relaxational density of states. — Before considering 
driving, we study the system's relaxational modes and 
rates. For a = 0, Eq. ^ is a generalized eigenvalue 
equation; its eigenvalues and -vectors are the relaxational 
rates and modes; Fig. [TJa gives an example. When de- 
formed along the mode |s„), the system relaxes expo- 
nentially to equilibrium with a rate given by the eigen- 
value s n < [T7j. We use s to refer to \s\ whenever no 
ambiguity results. 
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FIG. 2: Relaxational density of states. As the excess coordi- 
nation number Az (legend) vanishes, D(s) develops a power 
law divergence. The dashed curve ~ 1/s 0,47 . Inset: The 
low s data is collapsed by rescaling s with the crossover rate 
s* ~ Az 1 ' 9 . Exponents are best fits (see text). 



Fig. [2] displays the relaxational density of states D(s), 
the distribution of rates, averaged over approximately 
50 packings of N = 1024 disks for multiple values of 
Az. Several features stand out. First, close to jam- 
ming (Az — > 0), D(s) develops a power law divergence: 
packings at the transition possess many slowly relaxing 
modes. Second, the divergence is cut off at a crossover 
rate s* that vanishes with Az, so that, as packings are 
prepared closer and closer to the transition, the weight 
of slow modes grows. Finally, the fastest relaxations are 
on the order of the bare rate k/b = 1. 

Scaling of D(s). — To explain the characteristic fea- 
tures of D(s), we generalize the "cutting argument" of 
Wyart, Nagel, and Witten (WNW) to overdamped dy- 
namics. In seminal work, WNW introduced this varia- 
tional argument to explain the low frequency plateau in 
the vibrational density of states of undamped packings, 
D(Cl) |18| . Here we sketch relevant details. 

The WNW argument involves constructing a set of 
trial modes and estimating their frequencies, from which 
the scaling of D(fl) can be inferred. The trial modes 
are made by first "cutting" contacts at the boundary of 
patches of size t - this introduces floppy modes, zero fre- 
quency collective modes that involve no relative normal 
motions between disks - and then "stitching up" the cut 
with a sinusoidal envelope with wavenumber q ~ l/£. Be- 
cause of the sinusoid there are small relative normal mo- 
tions in the trial mode. Trial modes can be constructed 
for wavenumbers q > q* ~ Az, and their probability 
density is flat, D(q) ~ q° for q> q* . 

To generalize these results to damped dynamics, there 
are two key observations. First, we anticipate that, for 
sufficiently low frequencies, disk motions smoothly ap- 
proach their quasistatic values. Thus for long time scales 
the trial modes are also good trial modes for the over- 
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damped dynamics. Second, D(s) can be related to D(q) 
by inferring the "dispersion relation" , i.e. how s scales 
with q. If s ~ q 1 /", then D(s) ~ q°\dq/ds\ ~ s*" 1 . 

To estimate the overdamped dispersion relation, we as- 
sume that a trial mode |qwnw) is an approximate solu- 
tion to the homogeneous equation of motion, K\qwNw} + 
sB | <7wnw) ~ 0. Taking an inner product and bearing in 
mind Eqs. ([!]) and ^ for V and R, at the transition 
(z = z c ) one finds 
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where for scaling we refer to typical (relative) displace- 
ments. Normal motions in a trial mode come from the 
sinusoidal modulation: Am!| vnw ~ 9uwnw ~ g«wNW- 
In contrast, tangential motions are predominantly con- 
tributed by the floppy mode - they persist if q is sent 
to zero - so to leading order they are independent of q, 
Au 
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FIG. 3: Data collapse of the storage modulus G' (open 
symbols) and loss modulus G" (filled symbols) plotted for 
uo < 0.07 and varying distance to jamming Az (legend). The 



long (short) dashed curve indicates 



(w) scaling. Expo- 
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UWNW- In the limit of small q, then, the dis- nents are best fits (see text). Inset: G* prior to rescaling. 



persion relation is quadratic, s ~ q . Hence v = 1/2, and 
at jamming the density of states diverges as 
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(5) 



For z > z c , the dispersion relation remains quadratic for 

s > s* := (q*) 2 , or 
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so Eq. ^ continues to hold above s* . Eqs. Q and Q 
are our first main result. 

The predicted divergence and crossover scaling in D(s) 
can be tested by plotting s A D(s) versus s/Az x , as shown 
in the inset of Fig. [2j Eqs. ^ and ^ predict that the 
crossovers at low s will collapse for A = 2, and that above 
the crossover the curves will be flat for A = 1/2. We 
find the best collapse and plateau for A = 1.9(1) and 
A = 0.47(5), in good agreement with our predictions. 

The vanishing rate s* implies a diverging time scale 
t* := 1/s* ~ 1/ Az 2 . We stress that r* is different from 
the diverging time scale ~ 1/Az observed in undamped 
packings [18]. These two time scales come from differ- 
ent limiting cases of dynamics with both damping and 
inertia; when inertial effects dominate damping the dis- 
persion relation is linear rather than quadratic. 

Shear response. — Near jamming there is an abundance 
of slow relaxational modes. These must influence the re- 
sponse to shear (Fig. lip). We now show that the anoma- 
lous modes are responsible for shear thinning at frequen- 
cies to > s* , while for u < s* the response is quasistatic. 

For a sinusoidal strain with frequency u), the shear 
stress is a(t) = G* (u)j(t). Numerical results for G* on 
approach to jamming are plotted in the inset of Fig. [3] 
|19j . At high frequency cj > 1, loss dominates storage 
and both moduli reflect their microscopic counterparts: 
G" I to ~ G' ~ 0(1). This is because the packing cannot 



relax faster than the bare rate k/b and there is no inertia, 
hence high frequency deformations are affine and viscous 
stress dominates. Elasticity dominates for low frequen- 
cies, with the quasistatic modulus Go := lim w _>o G' di- 
minishing and the dynamic viscosity r) Q := lim^^o G" /lo 
growing on approach to jamming. Between these two ex- 
tremes is a range of frequencies over which G' and G" 
grow with sub-linear dependence on frequency, i.e. shear 
thinning; this range grows as jamming is approached. 

Scaling of G*{lo). — By expanding ^(s)) in the relax- 
ational modes and using the fact that the modes are B- 
orthogonal, (s m \B\s n ) — for m ^ n, Eq. ^ can be 
inverted to give an exact expression for G*: 
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with G*(u) = (|s„| + iw)(s„|i?|s„)/|(s n |7)| 2 . Such a re- 
ciprocal sum describes the modulus of a series circuit of 
viscoelastic elements. Further, the form of G* (w) is iden- 
tical to the shear modulus of a Kelvin- Voigt element with 
rate constant \s n \; each mode acts like a Kelvin- Voigt el- 
ement in the circuit (Fig. [T^i). We stress that Eq. ^ 
is not a model postulated ad hoc, but follows directly 
from inverting the equations of motion. It makes clear 
that linear rheology is controlled by the distribution of 
relaxation rates, D(s). 

The scaling of G* can be extracted from Eq. Q . Treat- 
ing each mode's shear component (s n |7) as a random 
variable independent of s„, near jamming one finds 
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Asymptotic analysis of the integral gives 

{Az ( uj/Az uj <s* 

uj 1 / 2 and G" ~ \ uj 1 ' 2 s* < u> < 1 (9) 
1 { UJ 1 < u. 

This is our second main result. 

Eq. Q provides (i) an analytical prediction for the 
scaling of the quasistatic shear modulus Go ~ Az, con- 
sistent with numerical and analytical modeling in the 
QS limit O |2Qj. It predicts (ii) a diverging viscosity 
t]q ~ 1/Az, which has not previously been predicted or 
measured. It shows (iii) that G' deviates qualitatively 
from the QS limit for uj > s* . Thus the QS approxima- 
tion is only reasonable for time scales much longer than 
t* ~ 1/Az 2 , which diverges at jamming. Finally, (iv), 
Eq. ([9]) predicts that the storage and loss moduli display a 
shear thinning regime scaling as uj a , with A = 1/2 being 
the same exponent that governs the divergence of D(s). 
The regime extends to frequencies as low as uj ~ s* , 
and hence to zero frequency at the jamming transition. 
Such scaling is therefore a critical effect, and the critical 
jammed state is a shear thinning complex fluid. 

Eq. ([9l also predicts that plotting the rescaled complex 
modulus G* / Az^ and rescaled frequency ujj Az x will pro- 
duce data collapse for fi = 1 and, again, A — 2. Indeed we 
find excellent collapse for fj, — 1.05(5) and A = 1.95(5) 
(Fig. [3j, in good agreement with the predicted values. 
The loss modulus is linear for low rescaled frequency 
(short dashed curve) , as predicted, while for larger values 
both G' and G" approach a power law scaling G* ~ («w) A 
(long dashed curve); the value A = fi/X = 0.54(4) again 
agrees well with its predicted value of 1/2. As Eq. ^ also 
captures the trivial high frequency scaling, its predictions 
are fully confirmed. While our numerics are restricted to 
two dimensions, the predicted exponents are independent 
of d, as are all known jamming exponents [3]. 

Outlook. — We have identified and explained the diver- 
gence and vanishing crossover in the relaxational density 
of states. We have also shown that slow relaxational 
modes determine the form of the complex shear modulus 
near jamming, including the quasistatic shear modulus 
Go and dynamic viscosity 770, the critical shear thinning 
regime, and the diverging relaxation time r* governing 
their crossover. 

The predicted QS scaling has yet to be observed ex- 
perimentally; low frequency rheology is complicated by 
aging processes such as coarsening, which are not mod- 
eled here [8]. By contrast, the anomalous shear thinning 
regime has already been observed [7J |5j. While prior 
work by Liu et al. also predicts w 1 / 2 scaling 7], it invokes 
the fluctuation-dissipation theorem and does not explain 
the connection to other scaling regimes or to jamming. 
Molecular dynamics simulations below jamming also find 
A ss 0.5, consistent with our prediction [12] . 

We have demonstrated that collective effects can give 
rise to anomalous uj 1 ! 2 scaling in a system where micro- 



scopic viscous forces are linear. One might also look for 
sources of nonlinearity in the microscopic interactions, 
but we know of no such mechanism common to liquid 
and organic foams, emulsions, and microgel suspensions, 
all of which display the same shear thinning rheology [5] . 
We therefore conclude that A = 1/2 is one of the few 
experimentally observed jamming exponents [5]. 
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